A generalized quantum microcanonical ensemble from random matrix product states 
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We propose a tensor network algorithm for the efficient sampling of quantum pure states belonging 
to a generalized microcanonical ensemble. The algorithm consists in an adaptation of the power 
method to a recently introduced ensemble of random matrix product states. The microcanonical 
ensemble that we consider is characterized by the fact that the participating energy eigenstates are 
not required to have identical statistical weight. To test the method we apply it to the Heisenberg 
model with an external magnetic field, and we find that the magnetization curves, due to the 
microcanonical constraint, are qualitatively different from those obtained in the canonical ensemble. 
Possible future applications include the study of isolated quantum systems evolving after a quantum 
quench. 
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Introduction. — The problem of representing statistical 
mechanics ensembles within the more fundamental quan- 
tum theory has recently attracted new interest both at 
the experimental and theoretical level [D-Q. On the ex- 
perimental side, the ability to isolate quantum systems 
from their surrounding triggers new questions on quan- 
tum thermalization and quantum equilibration dynamics 
[Tol | . On the theoretical side, a conceptually new ap- 
proach based on quantum typicality has emerged as an 
alternative paradigm for the understanding of quantum 
statistical mechanics and thermodynamics 0, H|. Since 
in quantum mechanics ensembles are represented by den- 
sity matrices, there are two ways in which they can be 
obtained: (i) after a partial trace over an environment 
which is entangled with the system in a global pure state, 
(ii) averaging over an ensemble of global pure states. 
The paradigmatic change provided by quantum typical- 
ity has taken advantage of the partial trace point of view. 
Loosely speaking quantum typicality exploits the concen- 
tration, around the average value, of some random vari- 
able in high dimensional spaces. This can be used for ex- 
ample to justify the occurrence of the canonical ensemble 
as the reduced density matrix of a typical random pure 
state 0, HI EH, E3 • On the other hand, averaging over en- 
sembles of microscopic configurations -point (ii) abovc- 
is yet another effective way of representing properties of 
physical systems [1, EH ■ I n this work we adopt this sec- 
ond approach, and we provide a tensor network algorithm 
to sample from a generalized quantum microcanonical en- 
semble, which differs from the standard microcanonical 
construction because the requirement for the participat- 
ing eigenstates to have equal statistical weights is relaxed 
to a more realistic one Q. The algorithm consists in 
an adaptation of the power method to Matrix Product 
States (MPS), and it exploits statistical properties of a 
recently introduced ensemble of Random MPS (RMPS) 
[3 E3 j which is the key feature allowing for the compu- 
tational efficiency of the method. In fact, any alternative 



procedure which considers Haar-distributcd random pure 
states would be far less efficient than the present scheme, 
due to the fact that those states arc computationally 
hard to generate. For the Heisenberg model with an ex- 
ternal magnetic field, the qualitative comparison of the 
results of our simulations with a different approach ex- 
ploiting Haar-distributed states, also supported by exact 
calculations corroborates the validity of our method. 
The introduced computational procedure, and the rela- 
tive generalized microcanonical ensemble, can also be of 
relevance in the context of current experimental setups 
where the quantum system is almost isolated from the 
environment [H-Q. 

The RMPS ensemble. — MPS are quantum pure states 
fully specified by a set {A ai , i = 1, . . . , N} of relatively 
small dimensional matrices 11711: 
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Following the standard notation for open-boundary MPS 
we set the row and column indices {ii, i./v} to 1; while all 
other indices belong to the set of integers from 1 to x> the 
latter being the so called bond-dimension of the MPS. For 
two levels systems, it follows that an MPS is specified by 
no more than 2Nx 2 numbers which, for x n ot too large, 
is exponentially smaller than the 2^ parameters required 
by a typical quantum state in the same Hilbcrt space T-L. 
Indeed, any state can be written as a MPS for x large 
enough, but the advantage in using MPSs occurs only 
when x ~ poly (AT). Moreover, MPS are physically rele- 
vant since they are good approximations of ground states 
of short-ranged gapped Hamiltonians, and they can also 
be used to construct thermal states HI 17, 18 1. With 
the aim to better understand the usefulness of typical- 
ity in quantum statistical mechanics, in 1A\ an ensem- 
ble of random MPS has been proposed sharing some of 
the statistical features of Haar-distributcd random pure 
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Figure 1. Gaussian fits for the histograms of the expecta- 
tion values of the transverse field Ising chain with respect to 
RMPSs. For ease of exposition, in the main figure only the 
normal fitting curves are shown, while the inset shows a par- 
ticular histogram with the relative Gaussian fit, for x — 8- 
N = 10, x = {4, 16, 64}, 1000 realizations. 



states. According to the construction presented in [lij ]. 
a random MPS is denned considering the ^-dimensional 
matrices A aie ^^ as contiguous sub-blocks of a set of 
N independent Haar-distributcd unitarics Ui, each uni- 
tary having dimension 2\ (where 2, in general, is the 
dimension d of the local Hilbert space). In particular, 
a sufficient condition for the occurrence of typicality in 
the case of RMPS has been derived in [lij], and the exact 
calculation of the average density matrix of the ensem- 
ble has been obtained in [l5(. Interestingly enough, the 
average random MPS coincides with the one obtained 
from Haar-distributed random states, i.e. it is the max- 
imally mixed state. For an alternative construction of 
random physical states, though not related to MPS, see 
also In this letter we show how the RMPS en- 

semble, together with efficient and relatively simple com- 
putational techniques, can be used to sample from a gen- 
eralized quantum microcanonical ensemble introduced by 
Reimann Computational efficiency is inherited from 
well-known MPS features, while on the statistical side 
we also exploit some of the properties characterizing the 
RMPS ensemble. In particular we will use the fact that 
the average state is equal to the completely mixed state 
[\ip) (V'Oave = V^), where \tp) is a random MPS as in EqQ] 
and D is the dimension of T~L. One of the implications of 
the above identity is numerically checked in FigQ], where 
the expectation value of the energy of the transverse field 
Ising chain is sampled over 1000 realizations of indepen- 
dent random MPS. The results of the simulations confirm 
that, irrespective of the auxiliary bond dimension \i the 
average energy is given by [{tp\H \tp)] avc = TrH/D = 0. 
Moreover in FigQ] one observes that the variance of the 
distribution depends on \- We are not able to derive 
an analytic expression for the second moment of the dis- 
tribution of RMPS, which would provide the variance in 



that the variance of the energy distribution decreases as 
a power-law with increasing \i consistently with typi- 
cality properties of the ensemble of RMPS [lij]. In the 
case of Haar-distributed random states the characteriza- 
tion of the probability distribution for general observ- 
ables can be found in [22 - 28 1 . We shall now present a 
computational procedure allowing us to modify the en- 
ergy distribution from the initial Gaussian centred in to 
another normal distribution with a desired average value, 
and whose variance can be reduced with a good degree 
of control. 

A generalized microcanonical ensemble from RMPS. — 
Our goal is to modify the initial distribution of Hamilto- 
nian expectation values in such a way that we can arbi- 
trarily choose its center E, and further decrease its vari- 
ance 7. To accomplish this we use an iterative technique 
which, starting from an initial RMPS \ip) , prepares a 
final MPS \ip) { whose population = | (Ei\ip)f | 2 - with 
respect to the eigenstates \Ei) of the Hamiltonian - is 
concentrated on the eigenvectors close in energy to a 
given value E. The iterative technique is provided by 
the power method: consider an operator G with max- 
imum eigenvalue A such that G\X) = A|A). It is easy 
to show that any initial state \tp} , which is not orthog- 
onal to I A), can be used as a starting point to obtain a 
good approximation to |A), simply iterating the following 
operation 
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For our purposes we set G = 
£ is a specified energy value, 



I - ( i± ^ A ) , where 
and a is a parame- 
ter which depends on the spectral range of the local 
Hamiltonian H, and which has to be chosen so that 
G is positive semi-definite (a simple estimate is pro- 
vided by o — 2N5 + \E\, where S is the greatest ab- 
solute value of the Hamiltonian parameters). The power 
method will bring an initial RMPS \ip) closer to the 
eigenstate \E), satisfying H\E) = E\E). Indeed, af- 
ter a number k of iterations \ip) becomes — 
G k \4>)o/\\G k |^)o Ik- In the following we provide a sim- 
ple argument showing that the asymptotic distribution 
of energies is approximated by a Gaussian centred in 
E. From the last equation we can write the aver- 
age density matrix at the fc-th iteration as follows: 
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G 2k , where we assumed that with high probability 
the denominator does not significantly depend on the 
initial state (as it occurs for Haar-distributed states 
Q), and we used the explicit expression for the av- 
erage state. In the limit of many iterations we can 
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the plots of Fig[TJ Nevertheless, we can check numerically is restricted to those eigenvalues sufficiently close to E. 
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From the above argument one can derive a polynomial 
upper-bound (in the size of the system) on the number k 
of iterations needed in order to have a sufficiently small 
variance (given by 4<r 2 /k): k ~ er 2 oc iV 2 . For any practi- 
cal purposes though, our numerical simulations show that 
already one hundred iterations are sufficient to obtain en- 
ergy distributions that are concentrated enough, and that 
can be used to analyse the system. Fig. [2] shows the his- 
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Figure 2. (Color online). Energy histograms and their nor- 
mal fitting curves at different steps of the power-method it- 
erations. N = 50, x = 32, E = -20, fOOO realizations. 

tograms of energy values obtained at different steps of the 
iteration. The simulation confirms that the distribution 
is indeed a Gaussian whose mean approaches asymptot- 
ically the value E, and whose variance decreases alge- 
braically with the number of iterations. The standard 
Dirac delta function representation of the microcanonical 
ensemble is approached in the limit of many iterations: 

[l^> (V'Oite ^ S( H ~ E )- Als0 ' n0tc that Ref - S 

provides a similar but more rigorous derivation of the 
average state in the context of Haar-distributed random 
states, showing that indeed the asymptotic final distri- 
bution of energies is a Gaussian centred in E. For small 
systems, using exact diagonalization we observe that at 
the level of a single RMPS realization the populations Pi 
of energy eigcnstates, which at the beginning is spread 
over a wide range of energy values, gets more and more 
concentrated around the chosen E value during the itera- 
tion process. This is consistent with the characterization 
of a generalized microcanonical ensemble, as proposed by 
Reimann in In particular, the standard requirement 
for the microcanonical ensemble, i.e. that any participat- 
ing eigenstate has the same probability in the ensemble, 
is weakened to more general non-uniform distributions of 
eigenstates in the small energy window, like those that 
we obtain. Reimann argues that these generalized mi- 
crocanonical energy distributions have to be considered 
more realistic realizations of the quantum microcanonical 
ensemble Q. Indeed, for a not too small system, energy 
levels are extremely close to each other and experimen- 



tally it would be extremely hard to prepare an equal su- 
perposition of energy eigcnstates. On the other hand it 
should be much simpler to require a distribution of the 
energy population with a sharp peak and a small variance 
(see also [llj and references therein). Note that, in princi- 
ple, for sufficiently many iterations, one single realization 
of a random MPS can provide an instance of the quantum 
microcanonical ensemble, along the lines of the construc- 
tion presented in [r| for Haar-distributed states. On the 
other hand, we can exploit the fact that sampling over 
man y R MPS improves the estimate of averaged quanti- 
ties [30 ] . and we do not need to wait for many iteration 
steps to terminate as long as we allow for a small finite 
width in the distribution of energies around their mean 
value. From this point of view the kind of generalized mi- 
crocanonical ensemble that we obtain is provided by an 
ensemble of random pure states which, when expressed in 
the energy eigenbasis of the Hamiltonian, have support 
on an energy window whose width can be controlled by 
the number of iterations, while the number of sampled 
states controls the accuracy in the estimation of aver- 
aged quantities. Moreover, the distribution obtained by 
sampling the Hamiltonian expectation values is not only 
meaningful from a statistical mechanics point of view in 
the regime of many iterations. In fact, the fluctuations in 
the energy value obtained at finite and not too large k can 
be interpreted as an energy exchange between the sam- 
pled system and a small finite-size environment, weekly 
interacting with the system. More about this interpre- 
tation can be found in j3f[ where the so-called Gaussian 
ensemble has been analysed in the context of classical 
statistical mechanics. For an alternative construction of 
a g eneralized quantum microcanonical ensemble see also 
32 . |33j , though in these works there is no constraint on 
the fact that, in the eigencnergies basis, the states in 
the ensemble have support only in a small energy win- 
dow, and in general they do not reproduce the standard 
canonical ensemble for subsystems. 

Application to spin chain Hamiltonians. — In this sec- 
tion we apply the RMPS power method to the Heisen- 
berg model with an external magnetic field: H = 

~ Eti I {of < +1 + ^f+i + of + ho\. To check 
the validity of the method we compare our results with 
similar quantities calculated for smaller systems in [lij ]. 
We compute the magnetization m z = -ZV -1 X)i=i K°f )Lve 
for the ferromagnetic model (J = 1), and the corre- 
lation function 4>{j) = N-^Zi [<°f CT f+i>] ave for the 
antifcrromagnetic model (J = —1). Both Fig. [3] and 
Fig. @] show lines labelled by different energy densities 
u = E/N, and each point on the curves is obtained 
averaging over 200 statistically independent RMPS (we 
checked that the standard deviations are always very 
small, of the order of 10 _3 J). We note that, from a 
computational point of view, in order to sample over 
many random MPS the algorithm can be trivially paral- 
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lclizcd to efficiently estimate statistical properties of the 
ensemble. The results of our simulations are qualitati vely 
and quantitatively very similar to those provided in , 
where the authors can also use data obtained from exact 
calculations in the canonical ensemble which are consis- 
tent with their findings. We note that, although our algo- 
rithm shares some aspects with the procedure presented 
in (Til, the two computational schemes are substantially 
different, mainly due to the fact that we make use of an 
MPS ensemble, allowing us to access much larger spin 
chains. From a physical point of view it is interesting to 
note that Fig. [3] shows a non-monotonic behavior of the 
magnetization with respect to the field, which is quite dif- 
ferent from what one obtains in the canonical ensemble at 
fixed temperature (i.e. monotonic magnetization curves) 
34|. Note that this difference between the two ensem- 
bles is evident when looking at the 'global' behavior of 
the magnetization curves, while 'locally' point-by-point 
one can always find canonical ensembles (characterized 
by different temperatures) providing the same results as 
the generalized microcanonical one [HI]. The difference 
is due to the energy constraint imposed by the quan- 
tum microcanonical ensemble, which does not allow the 
system to explore the same set of states as in the canon- 
ical case. In fact there are two main contributions to 
the Hamiltonian expectation value: one due to the in- 
teraction term, the other due to the magnetic field term. 
For small increasing values of the magnetic field the sys- 
tem can increase its magnetization, while keeping the 
energy constant, due to the effect of the interaction term 
which can balance the different energy contribution com- 
ing from the magnetic field term (this explains the left 
side in Fig. [3]). On the other hand, when the field is big 
enough, a further increase in the magnetization can not 
be compensated by the interaction term, implying that 
the microcanonical energy constraint would not be sat- 
isfied. The only states allowed in this regime are those 
for which the magnetization decreases with an increasing 
field. Considering state-of-art experiments with almost 
isolated quantum systems, it would be interesting to have 
an experimental confirmation of this behavior in those 
setup where it is possible to access the energy window of 
interest, to control the external field and to measure the 
magnetization. Fig. [4] provides information on the decay 
of the correlations and on the way in which they can be 
affected by an external magnetic field. Also in this case 
our simulations compare well with those provided in Rcf. 

Conclusions and future directions. — We have pre- 
sented a tensor network algorithm which can be used 
to sample states belonging to a given energy window. 
Statistical properties of the initial ensemble of random 
MPS, and the exploitation of tensor network structures 
make the algorithm suitable for the simulation of large 
one-dimensional systems. This algorithm allows us to 
test new ideas in the foundations of quantum statistical 
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Figure 3. (Color online). Magnetization curves of the ferro- 
magnetic Heisenberg chain with an external field, for different 
values of the energy density. N = 50, \ — 16. 
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Figure 4. (Color online). Average correlation function of the 
antiferromagnetic Heisenberg chain with an external field, for 
different values of the energy density. The figure on the left 
shows the results with h — 0, while the figure on the right 
shows the results for h — 0.4. N = 50, \ = 16. 



mechanics regarding a generalization of the quantum mi- 
crocanonical ensemble, as presented in With respect 
to this, in the future we plan to study in more details 
computational aspects connected to quantum typicality 
in our ensemble of random MPS. Our approach suggests 
also a new algorithm for the simulation of quantum sys- 
tems in the canonical ensemble, with some features sim- 
ilar to the algorithm proposed in 



131 ] . In the context 



of quantum computation it would also be interesting to 
investigate the computational complexity of a quantum 
algorithm able to sample from the microcanonical distri- 
bution. Along the lines of our construction, one should 
combine a subroutine implementing a quantum pseudo- 
random circuit [3~HI ] used to generate random quantum 
states, together with the iterative application of the op- 
erator G in EqJ5] using the technique developed in [36j . 
Alternatively one could think of adapting a similar quan- 
tum circuit into the framework of mixed states quantum 
computation [37i [HI . 
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Generation of random MPS 

The Hilbert space of a system of N qubits is expo- 
nentially big in N, hence typically it would require an 
exponential number of coefficients to characterize a state 
in a given basis. Matrix product states provide a so- 
lution to this problem, since they allow for a compact 
description of the states [l?], IHj]. MPS are represented 
as linear combinations of basis vectors whose coefficients 
are given in terms of products of matrices. For a system 
composed of N particles, with local Hilbert space dimen- 
sion d, an MPS state can be specified using only a number 
dN of A ai matrices (where i = 1, . . . , N), of dimension 
at most x ( a parameter known as the bond dimension). 
4153W u f^M :i 2') ^ N coefficients required for a typical state in 
Lhc Hilbert space are compressed into a set of at most 
dN\ 2 numbers, which makes computations much more 
efficient, but of course also restricts the class of possible 
states. Fortunately this restricted set can represent good 
approximations for many interesting states of relevance 
in condensed matter. More specifically, MPS are able to 
represent any state whose entanglement obeys the area 
law or only violates it logarithmic. 

With respect to the boundary conditions, there are 
two different constructions of MPS: with periodic or with 
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open boundaries. In the first case the state is written as 

(7 

while for open boundaries the first and last matrices are 
row and column vectors 



vectors. Therefore the ensemble of random MPS gener- 
ated in this way is appealing for different reasons: it has 
a natural operational meaning, the state is obtained in a 
canonical form, and it allows to use properties of the Haar 
distribution. Typicality and other properties of the en- 
semble of random MPS have been studied in 
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\A^...A^-^)\a 1 a 2 ...a N ). 



One can also consider the case where the A matrices are 
the same at each site, a so-called homogeneous MPS. 
Together with periodic boundary conditions this also 
gives a translational invariant state. Note that a general 
non- homogeneous MPS, with or without open bound- 
aries, can also describe a translational invariant state 
[39^ . There is no unique correspondence between the set 
of pure MPS states and the set of A matrices: differ- 
ent set of matrices may generate the same MPS. This 
gauge degree of freedom can be fixed using two canonical 
forms for the MPS: left-canonical MPS or right-canonical 
MPS. In the left-canonical form the A matrices satisfy 
the following condition: A aii A ai = I. In the right- 
canonical form the A matrices satisfy the following con- 
dition: A ai A' 7 ^ = I. When constructing an MPS in 
canonical form the boundary sites can violate the above 
conditions, but this simply means that the state is not 
normalized and it also provides an efficient way to eval- 
uate its norm (l3 |. 



In defining the ensemble of random MPS |14| we took 
inspiration from the sequential generation of MPS (39l . 
|40| . according to which a matrix product state can be 
seen as the outcome of a sequential interaction of unitary 
matrices between local physical degrees of freedom and 
an ancillary space. Along this line a natural way for 
defining a random MPS starts by considering a random 
Haar-distributed unitary matrix U filj of dimension d%, 
which can be seen as composed of blocks in the following 
way (for simplicity we now restrict to the case d = 2) 



U = 



M V 
N W 



where each block has dimension \- Identifying the block 
M and N with the two A ai=1 '° matrices associated to 
a local qubit immediately provides the building block 
for the construction of a left-canonical random MPS. On 
the other hand a right-canonical MPS can be obtained 
by identifying the two A ai=1 '° matrices with M and V. 
Similar considerations hold for the block column matri- 
ces V and W, or for the row column matrices N and 
W. A non-homogeneous random MPS can be iteratively 
constructed by sampling from a set of N independent 
random unitary matrices, each of dimension 2%. In the 
case of open boundary conditions the first and last A ma- 
trices are simply obtained from random Haar-distributed 



Computational cost of the algorithm 

Our algorithm is an instance of the simple and widely 
used power method, an iterative technique which does 
not require any diagonalization or singular value decom- 
position. Instead one has to deal only with matrix-vector 
multiplications, during which an operator G is repeatedly 
applied to an initial state \tp). 

In order to estimate the computational resources re- 
quired by our scheme, we start by defining the operations 
needed in its implementation: generation of a random 
matrix product state \ip); construction of the operator G; 
application of the operator G to \ip); if needed the eval- 
uation of some expectation values. Except for the gener- 
ation of random MPS, all other computational steps are 
already well discussed in the literature (see for example 
[13]), an d are summarized here just for completeness. 

The generation of the initial random MPS has been 
described in the previous section. To estimate the com- 
putational resources, in the worst case, we only need to 
multiply by N (the length of the chain) the cost of the 
generation of a random unitary of size d\, where d is the 
local Hilbcrt space dimension (2 in the case of qubits) 
and x is the maximum allowed bond-dimension of the 
MPS. Since x is a t most a polynomial in TV, the state 
can be generated with an amount of resources scaling 
polynomially in TV. 

The construction of the operator G = I — ( H Z E ) ~ 
where H is the Hamiltonian of the system, while E and a 
are two given numbers- involves the Matrix Product Op- 
erator (MPO) representation of the identity matrix and 
of the first and second powers of the Hamiltonian. The 
MPO representation is the most effective way of dealing 
with operators in the MPS framework: it simply consists 
in the generalization to operators of the matrix product 
decomposition for states (see for example [13 ] for more 
details) 



O 



W° 



.W N ^ \a) (<r'|, (3) 



where the bond-dimension of the W matrices is typically 
small (equal to 5 or 9 as we will see). The MPO repre- 
sentation of the identity matrix is very simple, while less 
trivial is the representation of the Hamiltonian. An effi- 
cient way of constructing it is explained in [13] (see page 
142). It requires an MPO of bond-dimension 5 for the 
representation of H , while the cost of the representation 
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of H 2 can be optimized using a bond-dimension equal to 
9. We note that clearly the construction of the operator 
G has to be done only once at the beginning of the algo- 
rithm. Hence the MPO representation of G requires the 
storage of a number AN of matrices of dimension at most 
9 (25 if one does not optimize the construction of H 2 ). 

The third step in the algorithm consists in the appli- 
cation of the MPO G to the MPS \tp) which, for an MPS 
of bond-dimension \ an d an MPO of bond-dimension 
Xw, costs 0(Nd 2 x 2 Xw)- Since G is composed by a lin- 
ear combination of different MPOs, and since after the 
application of G to the bond-dimension of the state 
increases by a factor given by the MPO bond-dimension, 



we need to compress the resulting state to one with a 
smaller bond-dimension. This can be done using a vari- 
ational technique costing at most 0(Nx 3 X\v) UM- 

The above operations can be repeated for the de- 
sired number of times during the execution of the power 
method algorithm, and one can keep track of the expec- 
tation values of relevant observables at a cost which is 
linear in N, at most cubic in x and at most quadratic in 
Xw [I3- 

To conclude we note that in order to sample many 
different initial random MPS the algorithm can be triv- 
ially parallelized without penalizing the efficiency of the 
scheme. 



